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Abstract — Norden E. Huang et al. had proposed and 
published the Hilbert-Huang Transform (HHT) concept 
correspondently in 1996, 1998 [1]. The HHT is a novel 
method for adaptive spectral analysis of non-linear and 
non-stationary signals. The HHT comprises two 
components: - the Huang Empirical Mode 

Decomposition (EMD), resulting in an adaptive data- 
derived basis of Intrinsic Mode functions (IMFs), and the 
Hilbert Spectral Analysis (HSA1) based on the Hilbert 
Transform for 1 -dimension (ID) applied to the EMD 
IMF’s outcome. Although paper [1] describes the HHT 
concept in great depth, it does not contain all needed 
methodology to implement the HHT computer code [2]. 
In 2004, Semion Kizhner and Karin Blank implemented 
the reference digital HHT real-time data processing 
system for ID (HHT-DPS Version 1.4) [3-4]. The case 
for 2-Dimension (2D) (HHT2) proved to be difficult due 
to the computational complexity of EMD for 2D (EMD2) 
and absence of a suitable Hilbert Transform for 2D 
spectral analysis (HSA2) [7], [9]. The real-time EMD2 
and HSA2 comprise the real-time HHT2. Kizhner 
completed the real-time EMD 2 and the HSA2 reference 
digital implementations respectively in 2013 & 2014 [5- 
6]. Still, the HHT2 outcome synthesis remains an active 
research area. This paper presents the initial concepts and 
preliminary results of HHT2-based synthesis and its 
application to processing of signals contaminated by 
Radio-Frequency Interference (RFI), as well as optical 
systems’ fringe detection and mitigation at design stage. 
The Soil Moisture Active Passive (SMAP mission 
(SMAP) carries a radiometer instrument that measures 
Earth soil moisture at LI frequency (1.4 GHz polarimetric 
- H, V, 3 rd and 4 th Stokes parameters). There is abundant 
RFI at LI and because soil moisture is a strategic 
parameter, it is important to be able to recover the RFI- 
contaminated measurement samples (15% of telemetry). 
State-of-the-art only allows RFI detection and removes 


RFI-contaminated measurements. The HHT-based 
analysis and synthesis facilitates recovery of 
measurements contaminated by all kinds of RFI, 
including jamming [7-8]. The fringes are inherent in 
optical systems and multi-layer complex contour 
expensive coatings are employed to remove the unwanted 
fringes. HHT2-based analysis allows test image 
decomposition to analyze and detect fringes, and HHT2- 
based synthesis of useful image. 
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1. Introduction 

The HHT generated considerable interest in the 
worldwide engineering and scientific community. 
Multitudes of both research and engineering papers were 
published reporting HHT interpretation and application 
results based on HHT algorithm [2], [3]. The HHT is the 
first direct method of spectral analysis of non-linear and 
non-stationary signals since Napoleonic times. 
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Firstly, the HHT EMD component allows decomposition 
of a non-linear and non-stationary signal s(t), under most 
general conditions, into a finite set of IMFs. The IMF is a 
generalization of heritage transforms’ decomposition 
artificial basis such as polynomial terms in Taylor 
expansion, sine and cosine trigonometric functions basis 
in Fourier Transform or wavelets in the Wavelet 
Transform. Unlike its century old and most recent 
predecessors, the HHT’s EMD produces an outcome 
adaptive IMFs basis from the input data itself without 
any a priori assumptions about the data. 

Secondly, the Hilbert Transform (H) is applied to each 
IMF (2) yielding the complex part of the analytic signal 
defined by Dennis Gabor in [7]. 

The HHT IMFs basis is obtained by the following 
procedure called EMD, and EMD1 for the ID input 
vector. 

(a) Finding the local maxima and local minima 
points for the input function s(t) 

(b) Predicting these (a) sets of points outside the 
input function domain boundaries 

(c) Building the input function upper envelope Uen 
out of the local maxima and its predicted points 

(d) Constructing lower envelope Len out of the local 
minima and its predicted points (like a, b, c) 

(e) Computing the Median = (Uen + Len)/2 

(f) Finding the first IMF candidate by computing the 
difference IMFl=s(t) - Median over the input 
function domain. 

(g) Checking the IMF against the IMF criteria - 

(h) An IMF satisfies two conditions: 

- that the number of extrema and the number of 
zero crossings must either equal or differ at 

most by one in the whole data set and 

- that the mean value of the envelope defined by 
the local maximum and the envelope defined by 
the local minima at any point is zero 

(i) Repeating the procedure beginning at (a) by 
setting s(t) = IMF or sifting. 

Denoting a specific IMF function as d(t), the analytic 
signal and its spectral characteristics are defined by 
equations (l)-(6): analytic signal for d(t) (1) Hilbert 
Transform as convolution (*) of d(t) and l/7it (2), analytic 
signal representation in polar coordinates (3) allows to 
obtain instantaneous amplitude of the analytic signal (4), 
its local phase (5) and instantaneous frequency (6). 
Equations (1) - (6) form the HHT theoretical foundation 
for HHT1: 

a(t) = d(t) + iH(d(t)) (1) 
d H (t) = d(t) * 1/ TCt (2) 

a(t) = A(t)e i 0(t) 


A(a(t)) = V d 2 + d 2 H 

(4) 

0(t) = arctan(dH / d) 

(5) 

f(t) = (1/271 )0(t)’ 

(6) 


The HHT was selected as the NASA invention of the year 
2003 and its digital implementation, the HHT-DPS VI. 4, 
received the Regional Excellence in Technology Transfer 
Award in 2006 from the Federal Laboratory Consortium 
Mid-Atlantic Region. The HHT-DPS implemented the 
HHT for ID and became a real-time reference system 
with which all other implementations and results could be 
compared. It was maintained at the Goddard Space Flight 
Center (GSFC) where it was commercialized and licensed 
to DynaDx Corporation in 2009, leading to its 
development into a commercial product. Within the 
GSFC Internal Research and Development Program 
(IRAD), the HHT was under continuous research which 
resulted in development of the real-time EMD2 and 
HSA2 prototype technologies in 2013 and 2014 
correspondently. The complexity of HHT2 development 
required abstraction from applications to some degree. 
Presently, in view of the HHT2 completion, we are 
turning our attention to the HHT2 applications, namely 
the synthesis based on HHT2. Syneren Technologies 
Corporation further implemented the EMD2 as a 
reference real-time system commercial product and its 
applications are being explored under the 2013-2014 
Space Act Agreement between NASA/GSFC and Syneren 
Technologies Corporation. These applications include 
hyper-spectral image analysis, fringe detection and 
mitigation and RFI noise detection and science data 
recovery using HHT2-based synthesis. 

Following in Figures 1-2 is the presentation of HHT1 
EMD1 application to the analysis of length-of-day file 
LOD78.csv comprising 3250 data points that were 
collected in 1978 by the Jet Propulsion Laboratory (JPL). 



Figure 1. Graph of input file LOD78.csv 



(3) 
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Figure 2. HHT1 EMD1 IMFs graph for length- 


of-day variations in microseconds LOD78.csv 

The 7 IMFs and 1 residual were obtained running HHT- 
DPS version 1.4 on a Lenovo PC under Windows 7 OS 
and EMD1 configuration with such parameters: Sifting 
Criteria 3, Maximum IMF 100, Maximum # Sifts 50,000, 
Sift Type - Extrema Sift, Endpoint Selection - Pattern 
Prediction. 

2. Reference Real-Time EMD2 

The HHT2, similar to HHT1, comprises two modules the 

EMD2 and HSA2: 

HHT2 = {EMD2, HSA2} 

The problem of HHT2 digital implementation is two-fold: 

- EMD2 computational complexity has a formidable 
order of 0(N 4 ), which for an NxN image, where N=2048, 
is O(10 14 ). It is complex because it involves constructing 
image surface envelopes, envelope prediction at 4(N - 1) 
boundary points, and envelope smoothing. 

- HSA2 requires a suitable Hilbert Transform for 2D, 
which doesn’t exist in a form that allows reconstructing 
the 2D analytic signal. 

The computational complexity of EMD2 was resolved by 
developing its algorithm in a way that allows parallelism 
in execution. The input NxN image is then broken down 
into 64 smaller sub-arrays (Figure 2) and run on an 
NVIDIA GTX 690 Card. The EMD2 run time for a 2K x 
2K image with the number of Bidimensional Intrinsic 
Mode Functions (BIMF) set to 10, is a fraction of a 
second. The run time also depends on the number of 
extrema points in the input image. The performance 
evaluation was observed using image of raindrops on a 
lake (Figure 3), and the first BIMF (BIMF1) is depicted in 
Figure 4. 



Figure 2. EMD2 parallelism implemented 
on the NVIDIA GTX 690 PC/Win7 Card 

The NVIDIA card code was developed in C++ and the 
640 x 640 images in Figure 4 were obtained by Syneren 
Technologies Corporation with a run time of 0.1 seconds. 



Figure 3. Lake raindrops. Figure 4. The BIMF1. 


3. Reference Real-Time HSA2 

Spacecraft advanced passive microwave radiometry 
makes use of natural thermal emissions to remotely sense 
Earth phenomena of interest to science in the L-band 
(1.4GHz surface moisture, for example). Here, a 
terrestrial signal-to-space suffers less attenuation by the 
intervening canopy or atmosphere. Unfortunately, relative 
insensitivity of the L-band region to atmospheric effects 
also makes it an extremely attractive spectral range for 
wireless radar & cell communications that are causing 
RFI with the spaceflight science radiometer instrument’s 
terrestrial phenomenon signal of interest. State-of-the-art 
technology detects an RFI event and deletes the entire 
contaminated measurement point. Excision of just the 
RFI component from the composite signal is presently a 
tall challenge. The three goals described in this paper 
related to research and development of the HHT2 HSA2 
components are: 

1) Develop the engineering tool EMD2 for a direct 
decomposition of the RFI-contaminated signal into 2D 
(Bi-dimensional) IMFs (BIMFs) in real-time and use 
IMFs and BIMFs to recover RFI-contaminated science 
data points 

2) Research and develop the Hilbert-Huang Transform for 
2D (HHT2) last component - the Hilbert Spectral 
Analysis for 2D (HSA2) and complete the HHT2 theory. 
The HSA2 is required for reconstructing the 2D analytic 
signal (monogenic signal or monogenic) from which the 
2D image instantaneous amplitude and phase are 
obtained. 

3) Verify the new technology tool prototype by a 
controlled drone experiment at Wallops Flight Facility, 
comprising a miniature radiometer, a controlled RFI 
source, and science data processing on the ground. 

We constructed a bi-dimensional By-orthant Hilbert 
Transform (BHT or Rietz Transform) digital 
implementation that allows construction of the 2D 
analytic signal a(x, t) for the input image f(x, y). This was 
accomplished by using a combination of directed Hilbert 
Transforms along the OX and OY axis or Hx and Hy, as 
well as the total Hilbert Transform Ht =HxHy of f(x, y): 

BHT(f(x, y)) = Hx + Hy(f(x, y) - Ht(Ht)) (7) 

Then 

a(x, t) = f(x, y) + i*BHT(f(x, y)) (8) 

The instantaneous amplitude, phase, and frequency are 
derived for 2D as follows: 

a(x, t) = A(x, t)e' 0(x, t) (9) 
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A(a(x, t)) = V d 2 -2dxdH+ d 2 h (10) 

0(x, t) = arctan(dH / d) (11) 

f(v) = (1/271 )0(x,t,v) (12) 

This completes the research and development of the 
HHT2, where: 

HHT2 = {EMD2, HSA2} (13) 

4. Heritage Decompositions and Analysis 

4.1 Taylor Series Functional Analysis 

In mathematics, a Taylor series is a representation of a 
function as an infinite sum of terms that are calculated 
from the values of the function's derivatives at a single 
point. A one-dimensional Taylor series is an expansion of 
a real function, f(x), about a point, x=a, is given by: 

f(x) = f(a) + f(a)(x-a) + f ’(a)(x-a) 2 /2! + ... + F n) (x-a)7n! 
+... (14) 

The analysis of f(x) Taylor series and synthesis of its 
approximation to f(x) is done in the polynomial terms 
domain and is straight forward. 

4.2 FFT1 and FFT2 Spectral Analysis 

The digital signal s(t) fftl outcome-based analysis is 
depicted by the following sequence: 

s(t) -> fftl(s(t)) -> frequency spectrum vector -> spectrum 
vector processing (filtering) -> g(t) (15) 

4.3 HHT1 and HHT2 Analysis 

The digital signal s(t) HHT1 outcome-based analysis is 
depicted by the following sequence: 

s(t) -> (16) 

EMD1 {IMFi(t), IMF 2 (t), ..., IMF k (t)} -> 

HSA1 {a(IMFi(t)), a(IMF 2 (t)), ..., a(IMF k (t))} -> 
Ai(h), A 2 (h), ..., A k (h) -> 
f!(ti),f2(ti),...,fk(ti) 

where fi(t) is the instantaneous spectrogram for IMFi(t) 
and Ai(t) is the instantaneous amplitude of the 
instantaneous frequency component fi(t) for l<=i<=k. 

The analysis, performed in HHT1 frequency domain, is 
similar to the analysis in frequency domain performed for 
heritage transforms, such as the fftl depicted above. 

It is important to notice that HHT1 instantaneous 
characteristics of a signal are all bounded to the time U 
vertical slice and its synthesis is directly based on this 
time slice h. 


5. Heritage FFT Based Synthesis 

The fft based synthesis comprises processing of the 
frequency spectrum vector resulting in frequency 
synthesized vector g(t). The synthesis of the 
approximation to the input signal s(t) is then obtained by 
applying the inverse Fourier Transform or ifft(g(t)). 

5.1 FFT1 Outcome Synthesis 

Given a time variable function s(t), the FFT1 is used over 
a sequence of short time periods, At, to construct two- 
dimensional spectro-temporal “scene”. This image is then 
processed using image-based filtering and the 
approximation to s(t) is synthesized using IFFT1. 

The case for 2D is described in the following Section 5.2. 

5.2 FFT2 Outcome Synthesis 

An example of an fft2 -based synthesis is presented in [9]. 
The concept of the two-dimensional spectro-temporal 
modulation filtering of the auditory model is implemented 
for the fft spectrogram (2-D signal analyses). A 2-D 
invert fft (fft2) is then used to synthesize the denoised 
audio signal. 

6. HHT2 EMD2 and HSA2 outcome Synthesis 
The HHT based synthesis is anchored in the fact that all 
spectral characteristics of the input signal are tied to the 
unique time slice t,. and after processing of the 
instantaneous spectrum, the synthesized approximation to 
the original input signal is obtained by direct selection of 
corresponding to h instantaneous frequency components 
such as Ai(h) for l<=l<=k, l<=i<n. 

The HHT based synthesis is two-fold: 

• EMD IMF qualitative level and 

• HSA spectrum quantitative level. 

The HHT-based synthesis concepts and terminology are 
as follows: 

• The HHT-based synthesis is time-slice U oriented 

• The HHT-based synthesis is direct and 
accomplished by 

• moving up from frequency component of interest 
gi(h) along the time-slice U to the analytic signal 

• ai(h)) amplitude and 

• selecting the IMFi(t) the IMFi(h) datum from 
which this frequency entity of interest originated 
from. 

Once these HHT-based synthesis’ concepts and the 
terminology are established, we can move to the next step 


of HHT applications, such as fringe phenomena analysis 
based on HHT2 analysis, and mitigation based on HHT2 
synthesis. 


6.1 EMD2 Outcome Based Synthesis 

The EMD based synthesis is performed in the IMFs time 
domain similar to that of Taylor series extension. The 
IMF components of qualitative interest are selected to 
represent components of the synthesized signal 
approximation 

s(t) -> (17) 

EMD1 {IMFi(t), IMF 2 (t), ..., IMF k (t)} -> 

{IMFu(t), IMFli 2 (t), ..., IMFi m (t) } 

6.2 HSA2 Spectrum Based Synthesis 

The development of the real-time EMD2 and the HSA2 
required a great deal of abstraction from the pressures of 
a specific application. With the completion of the real- 
time HHT2 implementation, it is time to turn to 
applications of HHT2 to real problems of spaceflight 
remote sensing data analysis and synthesis. 

7. HHT Based Applications 

Our interest is in applications related to remote sensing 
data processing and recognition of its man-made signal 
components, such as RFI signals. There are a few 
emerging methodologies in discriminating man-made 
signals - information theory based negative entropy, 
statistical kurtosis. Applications of HHT1 and HHT2- 
based analysis and synthesis to RFT-contaminated signal 
reconstruction is an important application for advanced 
radiometer remote sensing missions. The outline of the 
HHT-based synthesis for such an application is presented 
below in Section 7.1 and 7.2. 

7.1 HHT2 ANALYSIS OF RFI CONTAMINATED SIGNAL 
The HHT2-based RFI detection and mitigation is as 
follows: 

(1) Perform SMAP RFI Detection Algorithms on a SMAP 
measurements file fl and on-board statistically derived 
from oversampling products of size N (fmll, fml2, 
fml3, fml4) generating a vector RFI detect of size N with 
Os and Is, where indices to Is point to RFI contaminated 
samples. This is state-of-the-art in RFI detection and RF_ 
samples removal. Perform the HHT2 EMD1/EMD2 on 
SMAP measurements and products files of size N 
(fl, fmll, fml2, fml3, fml4} applying SMAP RFI 
algorithms to IMFs sets: 

{ IMFi(f l),IMFi(fml 1 ),IMFi(fml 2),IMFi(fm 13), 
IMFi(fm4)} for i=l to k. 

Synthesize the i-indices into a single RFI-index file and 
compare with above state-of-the-art. 


(2) RFI event persists in a few IMFs across coincidental 
maxima spatial locations. Determine, for each IMF 
vertical cross-cuts, the IMFs with coincidental maxima 
points. Combine them and subtract the total from input 
signal to recover science measurements. 

(3) Apply HSA2 to each BIMF of drone imaging 
radiometer and synthesize the RFI spectra (Figure 5). 

(4) Apply the Principal Component Analysis (PCA) to the 
k x n IMFs to synthesize the SMAP measurements 
contaminated by RFI. 

7.2 HHT2-based synthesis of RFI decontaminated 
signal 

The recovery of RFT-contaminated data points (synthesis) 
is accomplished by applying steps (l)-(4) above and the 
results are depicted below in Figures 5 and 6. 
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Figure 5. RFI events detected by HHT 

.. . . i m. „.i i .i iiii.i l 

o O O G> 

Figure 6. RFI spectrum with RFI-e vents removed 
and HHT spectrum with RFI-events recovered data. 

8. Synthesis of fringe decontaminated images 

In physics, interference is a phenomenon in which two 
waves superimpose to form a complex wave, and fringes 
are observed wherever the two or more waves overlap. 
The fringe analysis has a wide range of applications in 
physical and engineering measurement based on fringes 
by design (interferometers) and analyses of unwanted 
fringes aims at their detection and mitigation. Fringe 
analysis’ state-of-the-art processing is using linearity- 
based Fourier Transform that is inadequate due to fringe 
phenomenon non-linear and non-stationary nature. The 
breakthrough technology for fringe analysis presented in 
this paper is based on the HHT2, allowing fringe 
detection and mitigation and phase untangling in real- 
time by firmware assisted HHT2. 
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Figure 7. EMD2 input Figure 8. Fringe BIMF1 
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10, 2004, Big Sky Montana, USA 


We have described the HHT2 reference digital 
implementation in real-time and initial solution 
methodology for HHT2 -based monogenic signal analysis 
synthesis. The HHT2 is a reference real-time data 
processing system. It processes a 2K x 2K image under a 
second for number of BIMFs larger than 10. It is a 
reference system in respect that it provides the HHT2 
computer code and its hardware platform specification so 
that an application can be run on an equivalent system 
and the result compared to this reference system 
maintained by Syneren Technologies Corporation under a 
Space Act Agreement with NASA. We have developed 
the EMD2 algorithm that allows parallelism and 
implemented HSA2 the Rietz Transform. We applied the 
HH2-based analysis to RFI detection and RFI- 
contaminated science data synthesis. We demonstrated 
the HHT2 applicability to fringe patterns decomposition 
analysis and image synthesis. 
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